% This program calculates the moments in the model we attempt to match to
% the data
function dataset = generate_summary_statistics(Prob, Ss, rs, lambdas, tol, max_It);
    global bet0 bet1 a b w e pi_theta N_a N_w N_e w_max w_min alpha...
        pi_e Pr skill_shock r_base xi1 myiter rho var_phi...
        shares w_bar ss_calc unemp zeta myslope elasts nu xi eta mu_a thetas tau sigmae share log_e_sigma sig idx_pe

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Steady-state choices and distribution
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Pre-allocate and compute variables
    dataset = cell(5, 23); % Dataset for the various moments
    ln_e_hat = -1/2*log_e_sigma^2; % mean of the non-normalized idiosyncratic wage shock
    N_n = length(Ss); 
    w_bar = dot(w, sum(Prob,2));
    
    % Calculate choice probabilities
    [P_S, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);

    % Initial Distribution
    myDist0 = P_S .* Prob;
    myDist0  = sum(bsxfun(@times, myDist0 , reshape(pi_theta, [1,1,1,2])),4);
    
    F_es0 = F_es;
    w0 = w;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Invert the housing equation to get lambdas
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    lambdas =  inv_hs(shares, rs);
    % For numerical convenience, this is equal to (\lambda_s)^(1/\phi_s) in
    % the paper
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Calculate the wage cutoff associated with college
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    edu_vec = educ_mat(:);
    pi_joint_vec = myDist0(:);
    
    edu_share_col = 0.18; % use the 17.8% percentile edu level at S.S. as college
    [edu_dec, index_edu] = sort(edu_vec,'descend');
    pi_joint_dec = pi_joint_vec(index_edu); %probability of sorted education
    
    cdf_edu = cumsum(pi_joint_dec);

    [unique_edu_dec,index,~] = unique(edu_dec,'first','legacy');
    find_edu_college = @(x)interp1(unique_edu_dec,cdf_edu(index),x,[],'extrap')-edu_share_col;
    edu_college = fzero(find_edu_college,median(edu_dec)); % gives education level at the edu_share percentile 
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Dynamics
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    for t = 1:5
        myiter = t;
        
        % Generate shock to spillover
        if t == 2
            eta = skill_shock;    
            F_es0 = F_es0 * skill_shock; % Adjust parent's wage for shock
        end
        
        pi_w = sum(Prob, 2);
        cdf_w = cumsum(pi_w,1);
        w_bar = dot(w, pi_w);
        
        w0_1 = w0;
        w0 = w;
        
        w_max0 = w_max;
        w_min0 = w_min;       
        
        idx_pe = find(cdf_w>0.25,1);
        
        % Solve Dynamic equilibrium
        [Ss, rs] = find_spillovers_dyn(Ss, lambdas, rs, Prob, tol, max_It);
        [P_S, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs);
        N_size = squeeze(sum(bsxfun(@times, Prob .* P_S , reshape(pi_theta, [1,1,1,2])),[1,2,4]));



        % Adjust parent's wage for the skill shock but don't let it impact
        % their optimal decisions
        if t == 1
            Prob_kids = update_dist_dyn(Prob, Pr, Ss, rs, 1);
        else
            Prob_kids = update_dist_dyn(Prob, Pr, Ss, rs, 0);
        end
        
        % Calculate the marginal distribution of w
        pi_w = sum(Prob, 2);
        cdf_w = cumsum(pi_w,1);
        
        % Find the median
        find_w_median= @(x)interp1(w0,cdf_w,x,[],'extrap')-.5;
        w_median = fzero(find_w_median,(w_max+w_min)/2);
        
        b0 = b;
        
        w_m = dot(w,pi_w);
        w_v = dot(((w-w_m).^2),pi_w);
        sd_w = w_v^(1/2);
        cofv = sd_w / w_m;
        
        % Distribution by neighborhood
        myDist1 = P_S .* Prob;
        myDist1  = sum(bsxfun(@times, myDist1 , reshape(pi_theta, [1,1,1,2])),4);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Return to Eduction
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        share_C = sum((educ_mat >= edu_college) .* myDist1,1:2)./sum( myDist1,1:2);
        Elogw_C = (sum((educ_mat > edu_college) .* log(wprimeS) .* myDist1,1:3)/sum((educ_mat > edu_college) .* myDist1,1:3));
        Elogw_NC = (sum((educ_mat < edu_college) .* log(wprimeS) .* myDist1,1:3)/sum((educ_mat < edu_college) .* myDist1,1:3));
        if ~isnan(edu_college)
            total_educ_return = Elogw_C-Elogw_NC;
        else
            total_educ_return=nan;
        end
        
        col_shares = squeeze(sum((educ_mat >= edu_college) .* myDist1,1:2)./sum(myDist1,1:2))';
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Average Education and ability
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
        
        avg_ed = squeeze(sum(educ_mat.* myDist1,1:2)./sum( myDist1 ,1:2))';
        avg_edu_pop = squeeze(sum(educ_mat.* myDist1,1:3)./sum( myDist1 ,1:3))';
        for c=1:3
            wprimeS_avedu(:,:,c) = (unemp + (b + avg_edu_pop * eta * (a' .* (bet0 + bet1*Ss(c)).^xi)) .* f_w(w))  ;
        end

        avg_ability = squeeze(sum(a'.* myDist1,1:2)./sum( myDist1 ,1:2))';
        w_p25 = find(cdf_w>0.25,1);
        w_p50 = find(cdf_w>0.5,1);
        w_p75 = find(cdf_w>0.75,1);
        avg_ability_w25 = squeeze(sum(a'.* myDist1(w_p25,:,:),1:2)./sum( myDist1(w_p25,:,:) ,1:2))';
        avg_ability_w50 = squeeze(sum(a'.* myDist1(w_p50,:,:),1:2)./sum( myDist1(w_p50,:,:) ,1:2))';
        avg_ability_w75 = squeeze(sum(a'.* myDist1(w_p75,:,:),1:2)./sum( myDist1(w_p75,:,:) ,1:2))';

        avg_logability = squeeze(sum(log(a)'.* myDist1,1:2)./sum( myDist1 ,1:2))';

        a_mat = repmat(a',N_w,1,N_n);
        w_mat = repmat(w,1,N_a,N_n);
        
        avg_a_city = sum(a_mat.* myDist1,1:3);
        avg_w_city = sum(w_mat.* myDist1,1:3);
        prod_aw = sum(w_mat.*a_mat.* myDist1,1:3);
        
        avg_asq_city = sum((a_mat.^2).* myDist1,1:3);
        avg_wsq_city = sum((w_mat.^2).* myDist1,1:3);

        
        corr_wa = (prod_aw - avg_w_city*avg_a_city)/(sqrt(avg_asq_city-avg_a_city^2)*sqrt(avg_wsq_city-avg_w_city^2));

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Gini Index
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        my_prod = pi_w .* pi_w' .* abs(w0 - w0');
        gini = sum(my_prod(:))/(2*dot(pi_w, w0));
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Income Dissimilarity Index
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        rich_cutoff = 0.8;
        poor = rich_cutoff; 
        rich = 1 - rich_cutoff; 
        agg_cdf = cumsum(myDist1, 1);
        
        P_poor_S = zeros(length(Ss),1);
        P_rich_S = zeros(length(Ss),1);
        
        find_rich_w = @(x)interp1(w0,cdf_w,x,[],'extrap')-rich_cutoff;
        w_rich = fzero(find_rich_w,(w_max0+w_min0)/2);
        
        share_p = squeeze(sum((w0 < w_rich) .* myDist1,1:2))/ sum(squeeze(sum((w0 < w_rich) .* myDist1,1:2)));
        share_r = (squeeze(sum((w0 >= w_rich) .* myDist1,1:2)))/(sum(squeeze(sum((w0 >= w_rich) .* myDist1,1:2))));
        
        dissim = sum(abs(share_p - share_r))*0.5; 

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Ability Dissimilarity Index
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        ab_cutoff = 0.8;
        low_ab    = ab_cutoff; 
        high_ab = 1 - ab_cutoff; 
        agg_cdf = cumsum(myDist1, 2);
        pi_a = sum(Prob, 1);
        cdf_a = cumsum(pi_a',1);

        P_lowab_S = zeros(length(Ss),1);
        P_highab_S = zeros(length(Ss),1);
        
        find_ab = @(x)interp1(a,cdf_a,x,[],'extrap')-ab_cutoff;
        a_highab = fzero(find_ab,1);
        
        share_la = squeeze(sum((a < a_highab)' .* myDist1,1:2))/ sum(squeeze(sum((a < a_highab)' .* myDist1,1:2)));
        share_ha = (squeeze(sum((a >= a_highab)' .* myDist1,1:2)))/(sum(squeeze(sum((a >= a_highab)' .* myDist1,1:2))));
        
        dissim_ab = sum(abs(share_la - share_ha))*0.5; 
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Dissimilarity Index by Education
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        share_lowedu  = squeeze(sum((educ_mat < edu_college) .* myDist1,1:2))/ sum(squeeze(sum((educ_mat < edu_college) .* myDist1,1:2)));
        share_highedu = (squeeze(sum((educ_mat >= edu_college) .* myDist1,1:2)))/(sum(squeeze(sum((educ_mat >= edu_college) .* myDist1,1:2))));
        
        dissim_edu = sum(abs(share_lowedu - share_highedu))*0.5;
        
        rich_ed_avg = sum(educ_mat .* myDist1 .* (w > w_rich),1:4)/sum(myDist1 .* (w > w_rich),1:4);
        poor_ed_avg = sum(educ_mat .* myDist1 .* (w < w_rich),1:4)/sum(myDist1 .* (w < w_rich),1:4);
        wealth_ed_avg = [rich_ed_avg, poor_ed_avg];

      
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Rank-Rank
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% intergenerational income rank correlation
        N_bp1 = 125;
        q = linspace(0,1,N_bp1);
        w_q = zeros(1,N_bp1);
        w_q(1) = w_min0;
        w_q(N_bp1)= w_max0+1e-5;
        for i = 1:(N_bp1-2)
          find_q_w = @(x)interp1(w0,cdf_w,x,[],'extrap')-q(i+1);
          w_q(i+1) = fzero(find_q_w,(w_max0+w_min0)/2);
        end
        
        w_qk = zeros(1,N_bp1);
        w_qk(1) = w_min;
        w_qk(N_bp1)= w_max+1e-5;
        cdf_w_kids = cumsum(sum(Prob_kids,2));
        for i = 1:(N_bp1-2)
          find_q_w = @(x)interp1(w,cdf_w_kids,x,[],'extrap')-q(i+1);
          w_qk(i+1) = fzero(find_q_w,(w_max+w_min)/2);
        end
        
        parents_rank = zeros(N_w,N_a);
        for m = 1:(N_bp1-1)
            for i = 1:N_w
                if w_q(m)<=w0(i) && w_q(m+1)>w0(i)
                    parents_rank(i,:) = m;
                end
            end
        end

        wprime = unemp + (b0 + F_es) .* f_w(w0);
        wprime = repmat(wprime, [1 1 1 N_e 2]);
        tmp = repmat(e, [1 N_w N_a N_n 2]);
        tmp = permute(tmp, [2 3 4 1 5]);
       
        wprime = wprime .* tmp;
        
        [~, kids_rank] = histc(wprime,w_qk);
        pi_joint_w_a_e = repmat(myDist1,1,1,1,N_e) .* permute(repmat(pi_e, [1 N_w N_a N_n]), [2 3 4 1]);
        mu_prank = sum(sum(sum(sum(parents_rank.*pi_joint_w_a_e))));
        mu_krank = sum(sum(sum(sum(kids_rank.*pi_joint_w_a_e))));
        mu_krank = pi_theta*mu_krank(:);
        prank_std = sum(sum(sum(sum((parents_rank-mu_prank).^2.*pi_joint_w_a_e))));
        prank_std = sqrt(prank_std(:));
        krank_std = sum(sum(sum(sum((kids_rank-mu_krank).^2.*pi_joint_w_a_e))));
        krank_std = sqrt(pi_theta*krank_std(:));
        rank_cov = sum(sum(sum(sum((parents_rank - mu_prank).*(kids_rank-mu_krank).*pi_joint_w_a_e))));
        rank_cov = pi_theta*rank_cov(:);

        rank_corr = rank_cov/(prank_std*krank_std);
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Calculate the flow of movers
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s
        place_e_diff_mass = zeros(N_n,N_n);
        
        for c = 1:N_n
            Probm = zeros(N_w,N_a);
            for j = 1:N_a
                for i = 1:N_w
                    wage = w0_1(i);

                    % The wage shock is i.i.d. and lognormally distributed
                    % Take midpoints and allocate mass to closest wage grid
                    % point
                    t1 = unemp + (b   + squeeze(F_es0(i,j,c))')  .* f_w(wage);
                    midpts = log(w0) - log(t1);
                    pts = ([midpts(1,:);(midpts(2:end,:) + midpts(1:end-1,:))/2; Inf] - ln_e_hat)/log_e_sigma;
                    mean_prob = diff(normcdf(pts));
                    tmp = bsxfun(@times, squeeze(mean_prob), Pr(j,:));
                    pr1 = myDist0(i,j,c) * tmp;
                    Probm(:, :) = Probm(:, :) + pr1;

                end
            end
            % Where their kids currently live

            myDist1_sub = squeeze(sum(bsxfun(@times, Probm .* P_S , reshape(pi_theta, [1,1,1,2])),4));
     
            % Difference in expected wage for their children
            
            place_e_diff_mass(c,:) = sum(myDist1_sub, 1:2);  
        end
        share_of_movers = place_e_diff_mass;


        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Chetty-Hendren Spillover returns
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Find 25th and 75th percentile
        p_low = 0.25;
        [~,find_w_p25]= min(abs(cdf_w-p_low));
        w_p25 = w0(find_w_p25);
        
        [~,find_w_p75]= min(abs(cdf_w-0.75));
        w_p75 = w0(find_w_p75);
        
        % Average wage of child with parent at 25th percentile
        % Interpolate between grid points above and below percentile
        idx = zeros(length(w0),1);
        idx(find_w_p25) = 1;
       
        p_blm_1 = sum((idx .* myDist1 .* wprimeS), 1:3) ./ sum((idx .* myDist1), 1:3);
        if cdf_w(find_w_p25) < p_low
            idx = zeros(length(w0),1);
            idx(find_w_p25 + 1) = 1;
            p_blm_2 = sum((idx .* myDist1 .* wprimeS), 1:3) ./ sum((idx .* myDist1), 1:3);
            wfirst = (p_low - cdf_w(find_w_p25))/(cdf_w(find_w_p25+1) -  cdf_w(find_w_p25));
            wlast = (cdf_w(find_w_p25+1) - p_low)/(cdf_w(find_w_p25+1) -  cdf_w(find_w_p25));
        else
            idx = zeros(length(w0),1);
            idx(find_w_p25 - 1) = 1;
            p_blm_2 = sum((idx .* myDist1 .* wprimeS), 1:3) ./ sum((idx .* myDist1), 1:3);
            wfirst = (cdf_w(find_w_p25) - p_low)/(cdf_w(find_w_p25) -  cdf_w(find_w_p25-1));
            wlast = (p_low - cdf_w(find_w_p25-1))/(cdf_w(find_w_p25) -  cdf_w(find_w_p25-1));
        end
        p_blm = p_blm_1*wlast + p_blm_2 *wfirst;
        
        % Average wage of child with parent at 75th percentile
        % Interpolate between grid points above and below percentile
        idx = zeros(length(w0),1);
        idx(find_w_p75) = 1;
        p_abm_1 = sum((idx .* myDist1 .* wprimeS), 1:3) ./ sum((idx .* myDist1), 1:3);
        if cdf_w(find_w_p75) < 0.75
            idx = zeros(length(w0),1);
            idx(find_w_p75 + 1) = 1;
            p_abm_2 = sum((idx .* myDist1 .* wprimeS), 1:3) ./ sum((idx .* myDist1), 1:3);
            wfirst = (0.75 - cdf_w(find_w_p75))/(cdf_w(find_w_p75+1) -  cdf_w(find_w_p75));
            wlast = (cdf_w(find_w_p75+1) - 0.75)/(cdf_w(find_w_p75+1) -  cdf_w(find_w_p75));
        else
            idx = zeros(length(w0),1);
            idx(find_w_p75 - 1) = 1;
            p_abm_2 = sum((idx .* myDist1 .* wprimeS), 1:3) ./ sum((idx .* myDist1), 1:3);
            wfirst = (cdf_w(find_w_p75) - 0.75)/(cdf_w(find_w_p75) -  cdf_w(find_w_p75-1));
            wlast = (0.75 - cdf_w(find_w_p75-1))/(cdf_w(find_w_p75) -  cdf_w(find_w_p75-1));
        end
        p_abm = p_abm_1*wlast + p_abm_2*wfirst; 
        
        
        
        pctls = [w_p25; w_p75];
        pctls_idx = [find_w_p25; find_w_p75];
        target_pctl = [p_low; 0.75];
        avgw = [p_blm; p_abm];
        
        % Calculate the CH returns
        spillovers = zeros(2,1);
        spillovers_all = zeros(2,1);
        spillovers_avedu = zeros(2,1);
        spillovers_a = zeros(2,1);
        spillovers_eff = zeros(2,1);
        spillovers_eff_avedu = zeros(2,1);
        for iter = 1:length(pctls)
      
            place_e_diff = zeros(N_n,N_n);
            place_e_diff_a = zeros(N_n,N_n);
            place_e_diff_all = zeros(N_n,N_n);
            place_e_diff_edu = zeros(N_n,N_n);
            place_e_diff_mass = zeros(N_n,1);
            place_e_diff_mass_a = zeros(N_n,1);
            place_e_diff_mass_all = zeros(N_n,1);
            for c = 1:length(Ss)    
                % Previous generation's neighborhood
                X_Prob = myDist0; 
                
                Probm = zeros(N_w,N_a);
                counter_neigh = setdiff(1:length(Ss), c);
                Probm = zeros(N_w,N_a);
                for j = 1:N_a
                    for i = 1:N_w
                        wage = w0_1(i);

                        % The wage shock is i.i.d. and lognormally distributed
                        % Take midpoints and allocate mass to closest wage grid
                        % point
                        t1 = unemp + (b   + squeeze(F_es0(i,j,c))')  .* f_w(wage);
                        midpts = log(w0) - log(t1);
                        pts = ([midpts(1,:);(midpts(2:end,:) + midpts(1:end-1,:))/2; Inf] - ln_e_hat)/log_e_sigma;
                        mean_prob = diff(normcdf(pts));
                        tmp = bsxfun(@times, squeeze(mean_prob), Pr(j,:));
                        pr1 = myDist0(i,j,c) * tmp;
                        Probm(:, :) = Probm(:, :) + pr1;

                    end
                end
               
                midx = pctls_idx(iter);
                mtarget = target_pctl(iter);
    
                full_wage_diff = (wprimeS);
                
                % Take children of movers at midx percentile
                idx = zeros(length(w0),1);
                idx(midx) = 1;
                % Calculate their counterfactual wage in each neighborhood
                myDist1_sub_1 = squeeze(sum(bsxfun(@times, Probm .* P_S .* idx, reshape(pi_theta, [1,1,1,2])),4));
                mfirst = sum(full_wage_diff .* squeeze(sum(myDist1_sub_1(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_1(:,:, counter_neigh), 1:3);

                myDist1_sub_all = squeeze(sum(bsxfun(@times, Probm .* P_S , reshape(pi_theta, [1,1,1,2])),4));
                dest(c,:) = (squeeze(sum(myDist1_sub_1(:,:,counter_neigh),1:2)))';
                a_mat = repmat(a',N_w,1,N_n);
                w_mat = repmat(w,1,N_a,N_n);
                ab_dest(c,:) = (squeeze(sum(a_mat(:,:,counter_neigh).*myDist1_sub_1(:,:,counter_neigh),1:2)./sum(myDist1_sub_1(:,:,counter_neigh),1:2)))';
                edu_dest(c,:) = (squeeze(sum(educ_mat(:,:,counter_neigh).*myDist1_sub_1(:,:,counter_neigh),1:2)./sum(myDist1_sub_1(:,:,counter_neigh),1:2)))';
                inc_dest(c,:) = (squeeze(sum(w_mat(:,:,counter_neigh).*myDist1_sub_all(:,:,counter_neigh),1:2)./sum(myDist1_sub_all(:,:,counter_neigh),1:2)))';

                
                
                % Here we take the whole population (not just movers)
                mfirst_all = sum(full_wage_diff .* squeeze(sum(myDist1_sub_1(:,:, :),3)), 1:2) ./ sum(myDist1_sub_1(:,:, :), 1:3);
                
                % Here we consider the wage gain under average education
                mfirst_avedu = sum(wprimeS_avedu .* squeeze(sum(myDist1_sub_1(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_1(:,:, counter_neigh), 1:3);

                % Here we fix an ability level (25th perc, mean, 75perc)
                pi_a = sum(Prob,1)';
                cdf_a = cumsum(pi_a,1);
                a_bar = dot(a, pi_a);

                amidx = find(a>=a_bar,1); % use this for mean ability
                [~,amidx]= min(abs(cdf_a-0.25)); %here use  0.25 or 0.75

                idxa = zeros(1,length(a));
                idxa(amidx) = 1;
                myDist1_sub_1_a = squeeze(sum(bsxfun(@times, Probm .* P_S .* idx.*idxa, reshape(pi_theta, [1,1,1,2])),4));
                afirst = sum(full_wage_diff .* squeeze(sum(myDist1_sub_1_a(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_1_a(:,:, counter_neigh), 1:3);
                frac_movers(c,:) = squeeze(sum(myDist1_sub_1_a(:,:, :), 1:2)/sum(myDist1_sub_1_a(:,:, :), 1:3))';

                % Here we take only the origin neighborhood and the effective destination                
                for i_n=1:2
                    w_eff = squeeze(sum(full_wage_diff .* squeeze(sum(myDist1_sub_1(:,:, counter_neigh(i_n)),3)), 1:2) ./ sum(myDist1_sub_1(:,:, counter_neigh(i_n)), 1:3));
                    w_eff_1(i_n,:) = [w_eff(c) , w_eff(counter_neigh(i_n))];
                end
                             

                % origin is common to both groups
                if sum(myDist1_sub_1(:,:, counter_neigh(1)), 1:3)>0
                    w_eff_1_c  = (w_eff_1(1,1)*sum(myDist1_sub_1(:,:, counter_neigh(1)), 1:3) + w_eff_1(2,1)*sum(myDist1_sub_1(:,:, counter_neigh(2)), 1:3))./sum(myDist1_sub_1(:,:, counter_neigh), 1:3);
                else
                     w_eff_1_c  = w_eff_1(2,1);
                end
                mfirst_eff = [w_eff_1_c w_eff_1(1,2) w_eff_1(2,2)];

                % Here we take only the origin neighborhood and the effective destination + evaluate at average education             
                for i_n=1:2
                    w_eff_avedu = squeeze(sum(wprimeS_avedu .* squeeze(sum(myDist1_sub_1(:,:, counter_neigh(i_n)),3)), 1:2) ./ sum(myDist1_sub_1(:,:, counter_neigh(i_n)), 1:3));
                    w_eff_1_avedu(i_n,:) = [w_eff_avedu(c) , w_eff_avedu(counter_neigh(i_n))];
                end
                             

                % origin is common to both groups
                if sum(myDist1_sub_1(:,:, counter_neigh(1)), 1:3)>0
                    w_eff_1_c_avedu  = (w_eff_1_avedu(1,1)*sum(myDist1_sub_1(:,:, counter_neigh(1)), 1:3) + w_eff_1_avedu(2,1)*sum(myDist1_sub_1(:,:, counter_neigh(2)), 1:3))./sum(myDist1_sub_1(:,:, counter_neigh), 1:3);
                else
                     w_eff_1_c_avedu  = w_eff_1_avedu(2,1);
                end
                mfirst_eff_avedu = [w_eff_1_c_avedu w_eff_1_avedu(1,2) w_eff_1_avedu(2,2)];

                edu_first = sum(educ_mat .* squeeze(sum(myDist1_sub_1(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_1(:,:, counter_neigh), 1:3);
                ab_first = sum(repmat(a',N_w,1,N_n) .* squeeze(sum(myDist1_sub_1(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_1(:,:, counter_neigh), 1:3);

                % Again, interpolating:
                % calculate weights for nearest grid points
              
                if cdf_w(midx) < mtarget
                    idx = zeros(length(w0),1);
                    idx(midx + 1) = 1;
                    myDist1_sub_2 = squeeze(sum(bsxfun(@times, Probm .* P_S .* idx, reshape(pi_theta, [1,1,1,2])),4));
                    myDist1_sub_2_a = squeeze(sum(bsxfun(@times, Probm .* P_S .* idx.*idxa, reshape(pi_theta, [1,1,1,2])),4));
                    wfirst = (mtarget - cdf_w(midx))/(cdf_w(midx+1) -  cdf_w(midx));
                    wlast = (cdf_w(midx+1) - mtarget)/(cdf_w(midx+1) -  cdf_w(midx));
                else
                    idx = zeros(length(w0),1);
                    idx(midx - 1) = 1;
                    myDist1_sub_2 = squeeze(sum(bsxfun(@times, Probm .* P_S .* idx, reshape(pi_theta, [1,1,1,2])),4));
                    myDist1_sub_2_a = squeeze(sum(bsxfun(@times, Probm .* P_S .* idx .* idxa, reshape(pi_theta, [1,1,1,2])),4));
                    wfirst = (cdf_w(midx) - mtarget)/(cdf_w(midx) -  cdf_w(midx-1));
                    wlast = (mtarget - cdf_w(midx-1))/(cdf_w(midx) -  cdf_w(midx-1));
                end
               
                % Calculate value at other grid point
                mlast = sum(full_wage_diff .* squeeze(sum(myDist1_sub_2(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_2(:,:, counter_neigh), 1:3);
                mlast_all = sum(full_wage_diff .* squeeze(sum(myDist1_sub_2(:,:, :),3)), 1:2) ./ sum(myDist1_sub_2(:,:, :), 1:3);
                mlast_avedu = sum(wprimeS_avedu .* squeeze(sum(myDist1_sub_2(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_2(:,:, counter_neigh), 1:3);
                alast = sum(full_wage_diff .* squeeze(sum(myDist1_sub_2_a(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_2_a(:,:, counter_neigh), 1:3);
                edu_last = sum(educ_mat .* squeeze(sum(myDist1_sub_2(:,:, counter_neigh),3)), 1:2) ./ sum(myDist1_sub_2(:,:, counter_neigh), 1:3);
                for i_n=1:2
                    w_eff_l = squeeze(sum(full_wage_diff .* squeeze(sum(myDist1_sub_2(:,:, counter_neigh(i_n)),3)), 1:2) ./ sum(myDist1_sub_2(:,:, counter_neigh(i_n)), 1:3));
                    w_eff_1_l(i_n,:) = [w_eff_l(c) , w_eff_l(counter_neigh(i_n))];
                end
                
                if sum(myDist1_sub_2(:,:, counter_neigh(1)), 1:3)>0
                    w_eff_1_c_l  = (w_eff_1_l(1,1)*sum(myDist1_sub_2(:,:, counter_neigh(1)), 1:3) + w_eff_1_l(2,1)*sum(myDist1_sub_2(:,:, counter_neigh(2)), 1:3))./sum(myDist1_sub_2(:,:, counter_neigh), 1:3);
                else
                    w_eff_1_c_l  = w_eff_1_l(2,1);
                end

                mlast_eff = [w_eff_1_c_l w_eff_1_l(1,2) w_eff_1_l(2,2)];

                for i_n=1:2
                    w_eff_l_avedu = squeeze(sum(wprimeS_avedu .* squeeze(sum(myDist1_sub_2(:,:, counter_neigh(i_n)),3)), 1:2) ./ sum(myDist1_sub_2(:,:, counter_neigh(i_n)), 1:3));
                    w_eff_1_l_avedu(i_n,:) = [w_eff_l_avedu(c) , w_eff_l_avedu(counter_neigh(i_n))];
                end
                
                if sum(myDist1_sub_2(:,:, counter_neigh(1)), 1:3)>0
                    w_eff_1_c_l_avedu  = (w_eff_1_l_avedu(1,1)*sum(myDist1_sub_2(:,:, counter_neigh(1)), 1:3) + w_eff_1_l_avedu(2,1)*sum(myDist1_sub_2(:,:, counter_neigh(2)), 1:3))./sum(myDist1_sub_2(:,:, counter_neigh), 1:3);
                else
                    w_eff_1_c_l_avedu  = w_eff_1_l_avedu(2,1);
                end

                mlast_eff_avedu = [w_eff_1_c_l_avedu w_eff_1_l_avedu(1,2) w_eff_1_l_avedu(2,2)];
                
                % Take weighted average over grid points
                int_val = mfirst*wlast + mlast *wfirst;
                int_val_all = mfirst_all*wlast + mlast_all *wfirst;
                int_val_avedu = mfirst_avedu*wlast + mlast_avedu *wfirst;
                int_val_a = afirst*wlast + alast *wfirst;
                int_val_eff = mfirst_eff*wlast + mlast_eff *wfirst;
                int_val_eff_avedu = mfirst_eff_avedu*wlast + mlast_eff_avedu *wfirst;
                edu_val = edu_first*wlast + edu_last *wfirst;
                mass_int_val = wlast*myDist1_sub_1 + wfirst*myDist1_sub_2;
                mass_int_val_a = wlast*myDist1_sub_1_a + wfirst*myDist1_sub_2_a;
                
                % Average wage and wage mass
                place_e_diff(c, :) = int_val; 
                place_e_diff_all(c, :) = int_val_all; 
                place_e_diff_avedu(c, :) = int_val_avedu;
                place_e_diff_a(c, :) = int_val_a; 
                if c==1
                    place_e_diff_eff = int_val_eff; 
                else
                    place_e_diff_eff = [place_e_diff_eff int_val_eff];
                end
                if c==1
                    place_e_diff_eff_avedu = int_val_eff_avedu; 
                else
                    place_e_diff_eff_avedu = [place_e_diff_eff_avedu int_val_eff_avedu];
                end

                place_e_diff_edu(c, :) = edu_val; 
                place_e_diff_mass(c) = sum(mass_int_val(:,:, counter_neigh), 1:3);
                place_e_diff_mass_all(c) = sum(mass_int_val(:,:, :), 1:3);
                place_e_diff_mass_a(c) = sum(mass_int_val_a(:,:, counter_neigh), 1:3);
                place_e_diff(isnan(place_e_diff)) = 0;
                place_e_diff_all(isnan(place_e_diff_all)) = 0;
                place_e_diff_avedu(isnan(place_e_diff_avedu)) = 0;
                place_e_diff_a(isnan(place_e_diff_a)) = 0;
                
                place_e_diff_eff(isnan(place_e_diff_eff)) = [];
                place_e_diff_eff_avedu(isnan(place_e_diff_eff_avedu)) = [];
            end
            
            % Std deviation for average wage of moving
            std_dev = sqrt(var(sum(place_e_diff .* place_e_diff_mass,1)./sum(place_e_diff_mass,1)));
            std_dev_all = sqrt(var(sum(place_e_diff_all .* place_e_diff_mass_all,1)./sum(place_e_diff_mass_all,1)));
            std_dev_avedu = sqrt(var(sum(place_e_diff_avedu .* place_e_diff_mass,1)./sum(place_e_diff_mass,1)));
            std_dev_a = sqrt(var(sum(place_e_diff_a .* place_e_diff_mass_a,1)./sum(place_e_diff_mass,1)));
            std_dev_eff = sqrt(var(sum(place_e_diff_eff .* place_e_diff_mass,1)./sum(place_e_diff_mass,1)));
            std_dev_eff_avedu = sqrt(var(sum(place_e_diff_eff_avedu .* place_e_diff_mass,1)./sum(place_e_diff_mass,1)));
            
            if iter==1 
                place_e_diff_25 = place_e_diff;
                place_e_diff_all_25 = place_e_diff_all;
                place_e_diff_avedu_25 = place_e_diff_avedu;
                place_e_diff_a_25 = place_e_diff_a;
                place_e_diff_eff_25 = place_e_diff_eff;
                place_e_diff_edu_25 = place_e_diff_edu;
                place_e_diff_mass_25 = place_e_diff_mass;
                dest_25 = dest;
                ab_dest_25 = ab_dest;
                edu_dest_25 = edu_dest;
                inc_dest_25 = inc_dest;
                frac_movers_25 = frac_movers;
            else 
                place_e_diff_75 = place_e_diff;
                place_e_diff_all_75 = place_e_diff_all;
                place_e_diff_avedu_75 = place_e_diff_avedu;
                place_e_diff_a_75 = place_e_diff_a;
                place_e_diff_eff_75 = place_e_diff_eff;
                place_e_diff_edu_75 = place_e_diff_edu;
                place_e_diff_mass_75 = place_e_diff_mass;
                dest_75 = dest;
                ab_dest_75 = ab_dest;
                edu_dest_75 = edu_dest;
                inc_dest_75 = inc_dest;
                frac_movers_75 = frac_movers;
            end 

            % Spillover
            spillovers(iter) = std_dev/avgw(iter); 
            spillovers_all(iter) = std_dev_all/avgw(iter); 
            spillovers_avedu(iter) = std_dev_avedu/avgw(iter); 
            spillovers_a(iter) = std_dev_a/avgw(iter); 
            spillovers_eff(iter) = std_dev_eff/avgw(iter); 
            spillovers_eff_avedu(iter) = std_dev_eff_avedu/avgw(iter); 
            
        end


        

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% Intergenerational Mobility
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %% Upward Mobility
        % find  percentile (cdf_w was already defined above)
        qtils = [0.2, 0.4, 0.6, 0.8];
        ids_p = zeros(length(qtils),1);
        ids_c = zeros(length(qtils),1);
        for id = 1:length(qtils)
            [~,idx] = min(abs(cdf_w -qtils(id)));
            ids_p(id) = idx;
        end
        
        pi_w_kids = sum(Prob_kids,2);
        cdf_w_kids = cumsum(pi_w_kids);
        for id = 1:length(qtils)
            find_q_w_kids = @(x)interp1(w,cdf_w_kids,x,[],'extrap')-qtils(id);
            w_q_kids = fzero(find_q_w_kids,(w_max+w_min)/2);
            ids_c(id) = w_q_kids;
        end

        % find probability of staying below 25th percentile for each (ability,wage) pair
        Q1toQ = zeros(N_a,N_w, N_n);
        mobility = zeros(length(qtils),length(qtils));
        for ii = 1:length(qtils)+1
            for jj = 1:length(qtils)+1
                if ii == 1
                    mstart = 1;
                    mend = ids_p(1);
                elseif ii == length(qtils)+1
                    mstart = ids_p(length(qtils))+1;
                    mend = length(w);
                else
                    mstart = ids_p(ii-1)+1;
                    mend = ids_p(ii);
                end
                if jj == length(qtils)+1
                    w_c_low = ids_c(length(qtils));
                    w_c_high = max(w);
                elseif jj == 1
                    w_c_low = 0;
                    w_c_high = ids_c(1);
                else
                    w_c_low = ids_c(jj-1);
                    w_c_high = ids_c(jj);
                end
                Q1toQ = zeros(N_a,N_w, N_n);
                for t2=1:2
                    for i=1:N_a
                        for j=mstart:mend
                            for c = 1:N_n
                                my_share = Prob(j,i)*P_S(j,i,c,t2)*pi_theta(t2);
                                wj = w0(j);
                                fwj = (b0 + F_es(j,i,c)).*e.*f_w(wj);
                                    Q1toQ(i,j) = Q1toQ(i,j) + ...
                                        my_share*sum( pi_e((w_c_low < fwj ) & (fwj  <= w_c_high)));
                            end
                        end
                    end
                end
                if mstart ~= 1
                    mobility(ii,jj) = sum(Q1toQ(:))/(cdf_w(mend) - cdf_w(mstart-1));
                else
                    mobility(ii,jj) = sum(Q1toQ(:))/(cdf_w(mend) );
                end
            end
        end
        mobility;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Share Above and Below
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
         %% Upward Mobility
        % find 25th percentile (cdf_w was already defined above)
        [~,idx] = min(abs(cdf_w -.25));
        find_25_w = @(x)interp1(w0,cdf_w,x,[],'extrap')-0.25;
 
        pi_w_kids = sum(Prob_kids,2);
        cdf_w_kids = cumsum(pi_w_kids);
        find_25_w_kids = @(x)interp1(w,cdf_w_kids,x,[],'extrap')-0.25;
        w_25_kids = fzero(find_25_w_kids,(w_max+w_min)/2);

        % find probability of staying below 25th percentile for each (ability,wage) pair
        earn_more = zeros(N_a,N_w, N_n);
        
        for t2=1:2
            for i=1:N_a
                for j=1:N_w
                    for c = 1:N_n
                        my_share = Prob(j,i)*P_S(j,i,c,t2)*pi_theta(t2);
                        wj = w0(j);
                            earn_more(i,j) = earn_more(i,j) + ...
                                my_share*sum(pi_e((b0 + a(i)*F_es(j,i,c)).*e.*f_w(wj) >= wj));
                    end
                end
            end
        end
        earn_more = sum(earn_more(:));
          
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Income Ratio
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        tmp3 = squeeze(sum(myDist1 .* w0, 1:2)./sum(myDist1 ,1:2));
        avg_inc_ratio = tmp3;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Collect Results
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        dataset{t,1} = dissim;
        dataset{t,2} = gini;
        dataset{t,3} = avg_inc_ratio';
        dataset{t,4} = P_poor_S';
        dataset{t,6} = rank_corr;
        dataset{t,7} = N_size';
        dataset{t,8} = total_educ_return;
        dataset{t,9} = spillovers';
        dataset{t,10} = rs';
        dataset{t,11} = p_blm/p_abm;
        dataset{t,12} = cofv;
        dataset{t,13} = Ss';
        dataset{t,14} = (squeeze(sum((w0 >= w_rich) .* myDist1,1:2)))'  ./ N_size';
        dataset{t,15} = squeeze(share_C);
        dataset{t,16} = col_shares;
        dataset{t,17} = dissim_edu;
        dataset{t,18} =  sum((w0 < rs(3)) .* myDist1 ,1:4);
        dataset{t,19} = avg_ed;
        dataset{t,20} = avg_ability;
        dataset{t,21} = mobility;
        dataset{t,22} = share_of_movers;
        dataset{t,23} = earn_more;
        dataset{t,24} = wealth_ed_avg ;
        dataset{t,25} = place_e_diff_25;
        dataset{t,26} = place_e_diff_edu_25;
        dataset{t,27} = place_e_diff_mass_25;
        dataset{t,28} = place_e_diff_75;
        dataset{t,29} = place_e_diff_edu_75;
        dataset{t,30} = place_e_diff_mass_75;
        dataset{t,31} = dest_25;
        dataset{t,32} = ab_dest_25;
        dataset{t,33} = dest_75;
        dataset{t,34} = ab_dest_75;
        dataset{t,35} = inc_dest_25;
        dataset{t,36} = inc_dest_75;
        dataset{t,37} = spillovers_eff';
        dataset{t,38} = place_e_diff_eff_25;
        dataset{t,39} = place_e_diff_eff_75;
        dataset{t,40} = spillovers_a';
        dataset{t,41} = place_e_diff_a_25;
        dataset{t,42} = place_e_diff_a_75;
        dataset{t,43} = frac_movers_25;
        dataset{t,44} = frac_movers_75;
        dataset{t,45} = spillovers_all';
        dataset{t,46} = place_e_diff_all_25;
        dataset{t,47} = place_e_diff_all_75;
        dataset{t,48} = corr_wa;
        dataset{t,49} = dissim_ab;
        dataset{t,50} = edu_dest_25;
        dataset{t,51} = edu_dest_75;
        dataset{t,52} = avg_ability_w25;
        dataset{t,53} = avg_ability_w50;
        dataset{t,54} = avg_ability_w75;
        dataset{t,55} = spillovers_avedu';
        dataset{t,56} = spillovers_eff_avedu';
       
        myDist0 = myDist1;
        Prob = Prob_kids;
       
        F_es0 = F_es;
    end

end